function rest = indexJ(K,N,RNK)

%% generating all index 1--N, most stupid method
% if K>1
%     tmp=ones(N^(K-1),1);
%     oldleft =kron(tmp,indexJ(K-1,N));
%     newright=kron((1:N)',ones(size(oldleft,1)/N,1));
%     rest = [newright,oldleft];
% elseif K==1
%     rest= (1:N)';
% end

%% generating all index 0--N, most stupid method
% if K>1
%     tmp=ones((N+1)^(K-1),1);
%     oldleft =kron(tmp,indexJ(K-1,N+1));
%     newright=kron((0:N)',ones(size(oldleft,1)/(N+1),1));
%     rest = [newright,oldleft];
% elseif K==1
%     rest= (0:N)';
% end
%% generating all index 0--N, each sum is bounded by N
% if K>1
%     tmp=ones((N+1)^(K-1),1);
%     oldleft =kron(tmp,indexJ(K-1,N+1));
%     newright=kron((0:N)',ones(size(oldleft,1)/(N+1),1));
%     tmprest1 = [newright,oldleft];
%     tmprest2 = tmprest1(sum(tmprest1,2)<=N,:);
%     rest=unique(tmprest2,'rows');
%     
% elseif K==1
%     rest= (0:N)';
% end
%% generating all index 0--N, each sum is bounded by N, and control index r
r = adaptindex(RNK);
if K>1
    tmp=ones((N+1)^(K-1),1);
    oldleft =kron(tmp,indexJ(K-1,N,RNK));
    newright=kron((0:N)',ones(size(oldleft,1)/(N+1),1));
    tmprest1 = [newright,oldleft];
    tmprest2 = tmprest1(sum(tmprest1,2)<=N,:);
    clear temprest1
    tmprest3=unique(tmprest2,'rows');
    clear tmprest2
    rest=tmprest3(tmprest3(:,1)<=r(K),:);
    
elseif K==1
    tmprest= (0:N)';
    rest=tmprest(tmprest<=r(K));
    
end
